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Image Segmentation by Base Point Selection and Wavefront Propagation 
Field of the Invention 

[01] The present invention relates generally to images processing, and more 
particularly to partitioning an image into homogeneous segments. 

Background of the Invention 

[02] A large number of image segmentation methods are known. Those 
methods are generally region oriented, clustering oriented, or edge oriented. Those 
methods are either based on concepts of similarity or discontinuity of feature 
values, e.g., intensities or colors, of individual pixels (points) in the image. 

[03] A region linkage oriented segmentation method 'grows' regions that have 
homogeneous intensity or color features. That method starts from an initial 'seed' 
point or pixel. The region is extended by one pixel at a time. An adjacent pixel is 
added to the region when a feature distance between the adjacent pixel and the 
region is less than a predetermined threshold. The linkage method has usually no 
control over a shape of a current boundary of the region as determined by the 
thresholding. 

[04] Histogram oriented methods attempt to remedy a number of deficiencies 
of simple thresholding. Additionally, histogram segmentation compensates for 
shifts in mean intensity level because an image intensity histogram distribution is 
considered, rather than examining intensity values of individual pixels directly, as 
in region linkage above. However, histogram segmentation has no explicit notion 
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of connectivity. An implicit assumption is that pixels with similar intensities 
belong to the same region. This may not be true in general. 

[05] Edge-based segmentation takes into account that edges are usually 
boundaries between segments. Edges are detected by searching for local 
discontinuities in image features. Then, a tracing process is applied to link edges 
into continuous and connected segment boundaries. Edge linking techniques suffer 
from a number of problems. Not only must edge pixels be linked into edge lists, 
but these lists must also be linked so as to extract closed region boundaries. The 
closure problem is not straightforward to solve, because many edges may terminate 
in the same location or large gaps may need to be bridged between edges. Detected 
edges can often not form a set of closed curves, which surround connected regions. 

[06] It is desired to address the disadvantages of the above segmentation 
methods. 

[07] Level sets provide a general framework for modeling surface wavefront 
propagation, Sethian, "An analysis of flame propagation," Ph.D. Dissertation, 
Mathematics, University of Berkeley, 1982, and Sethian, "Curvature and 
evaluation of wavefronts," Commun. in Mathematical Physics, Vol. 101, pp. 487- 
499, 1985 

[08] Level sets have been used to track moving wavefront interfaces in a large 
number of different applications, from fluid mechanics to computer vision, Osher 
et al., "Fronts propagating with curvature dependent speed: Algorithms based on 
Hamilton-Jacobi formulation," Journal of Computational Physics, Vol. 79, pp. 12- 
49, 1988, Koepfler et al., "A multiscale algorithm for image segmentation by 
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variationar method," SIAM J. of Num. Analysis, 31-1, pp. 282-299, 1994, Malladi 
et al., "Shape modeling with wavefront propagation: a level set approach," IEEE 
Trans, on Pattern Analysis and Machine Intelligence, Vol. 17, pp. 158-175, 1995, 
Faugeras et al., "Variational Principles, Surface Evolution, PDE's, level set 
methods and the Stereo Problem," IEEE Transactions on Image Processing, Vol. 7, 
No. 3, pp 336-344, 1998, Siddiqi et al., "A Hamiltonian Approach to the Eikonal 
Equation," Workshop on Energy Minimization Methods in Computer Vision and 
Pattern Recognition, 1999, Droske et al., "A multilevel segmentation method," 
Proc. Vision, Modeling and Visualization, MPI Informatik, Germany, pp. 327-336, 
2000, Leventon et al., "Level Set Based Segmentation with Intensity and Curvature 
Priors," IEEE Workshop on Mathematical Methods in Biomedical Image Analysis, 
2000, Vemuri et al., "A Level-Set Based Approach to Image Registration," IEEE 
Workshop on Mathematical Methods in Biomedical Image Analysis, 2000, 
Paragios et al., "Geodesic active regions and level set methods for supervised 
texture segmentation" International Journal of Computer Vision, Vol.46, pp 223, 
2002, and Chan et al., "Variational PDE Models in Image Processing," Notices of 
the American Mathematical Society, 2003. 

[09] Prior to level sets, a motion of a curve, a boundary, or a 'wavefronf was 
represented typically by a discrete parameterization of an object as a set of points. 
Positions of the points are updated according to a given evolution equation using, 
for example, snakes, string methods, and marker particles. However, those 
Lagrangian techniques rely on a continual reparameterization of the surface as the 
surface becomes more complex. Those methods also tend to disintegrate as a 
topology of an evolving shape changes. 
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[010] Level sets use a totally different model to represent moving wavefronts. 
The level set method views the wavefront propagation as a Eulerian, initial value 
partial differential equation. The moving wavefront is regarded as a zero level set 
of higher dimensional function. The evaluation of this higher dimensional function 
resembles a Hamilton- Jacobi equation, which enables the forming of sharp surface 
gradients and the incorporation of curvature effects into the propagation of the 
wavefront. 

[Oil] Prior art image segmentation suffers from either requiring a priori 
information to initialize regions, being computationally complex, or failing to 
establish the color consistency and spatial connectivity at the same time. Therefore, 
it is desired to provide a method, based on level sets, for partitioning an image into 
homogeneous segments according to some consistency criterion, e.g., color or 
intensity. 

Summary of the Invention 

[012] The invention provides a method for partitioning an image into segments. 
Previous attempts at image segmentation either require a priori information to 
initialize segments, are computationally complex, or fail to establish concurrently a 
color consistency and spatial connectivity. Therefore, the method according to the 
invention uses a region growing method that is based on a level set. 

[013] The segmentation problem is represented as a monotonic wavefront 
propagation in an absorbing medium with varying speeds of the wavefront. 
Wavefronts are emitted iteratively from selected base points (pixels). At a selected 
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base point, a local variance of data reaches a minimum. This indicates that the base 
point is representative of a local neighborhood of the base point, 

[014] A hierarchical gradient operator is applied to determine a local variance. 
A speed at which the wave is propagated is determined by a color similarity of a 
point on a wavefront to a current coverage of the wavefront, and by edge 
information. Thus, the wavefront advances in an anisotropic spatial-color space. 
An absorbing function acts as a stopping criterion of the wavefront. 

[015] The invention takes advantage of a fast marching method to solve an 
Eikonal equation for finding travel times of the wavefronts. 

[016] The method according to the invention is superior to linkage-based 
region growing techniques because the method prevents 'leakage' and imposes 
compactness on the region without over-smoothing a boundary of the region. 

[017] Furthermore, the method according to the invention deals with sharp 
comers and changes in topology. The automatic segmentation method is Eulerian. 
Thus, the method is efficient. 

Brief Description of the Drawings 

[018] Figure 1 is a block diagram of a method for partitioning an image into 
segments according to the invention; 

[019] Figure 2 is a block diagram of further details of the method of Figure 1 ; 



5 



• a 



MERL-1562 
Porikli 



[020] Figure 3 is a block diagram of additional details of the method of Figure 

l; 

[021] Figure 4 shows an input image and a corresponding variance image; 

[022] Figure 5 shows an input image and a corresponding variance image; 

[023] Figure 6 shows an input image and a corresponding variance image; 

[024] Figure 7 is shows wavefront propagation for the image of Figure 5, and 
the corresponding boundary of the wavefront upon termination; 

[025] Figure 8 is an image partitioned into segments corresponding to the 
image of Figure 4; 

[026] Figure 9 is an image partitioned into segments corresponding to the 
image of Figure 5; and 

[027] Figure 10 is an image partitioned into segments corresponding to the 
image of Figure 6; 

Detailed Description of the Preferred Embodiment 
[028] Level Sets 

[029] A method for partitioning an image into homogeneous segments 
according to the invention uses a level set. The method combines a numerical 
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solution of hyperbolic conversation laws into a Hamilton-Jacobi setting. Given a 
moving closed hypersurface y/{ x, t) in which is a 'wavefronf in the 
segmentation process, an Eulerian formulation is produced for the motion of the 
hypersurface propagating along a direction of its normal with a speed F. The speed 
of propagation can be a function of various parameters, such as curvature, normal 
direction, position, etc. 

[030] The main idea behind the invention is to embed this propagating interface 
as a zero level set of a N +1 dimensional function 0. The function </>(x, 0) is 

defined by <p (x, 0) = ± k, where A: is a distance from a point x : (x, y) to a surface 

^(x, 0). Thus, an initial function </>(x, 0) : R" +l — > R has the property that 

^(x, 0) = [x | <p(x, 0) = 0], i.e., the hypersurface corresponds to points (pixels) that 

are at distances k = 0 from the surface or wavefront. 

[031] This means the zero level set of the evolving function <j) always matches 
the propagating hypersurface, i.e., (j> (x, t) = 0. By the chain rule, one can obtain 
(/), +V0(x,*)x' =0. (1) 

[032] Because an evolution of the wavefront, in its outward normal direction, is 
equal to the speed of the wavefront F = x'. | V(f> |,a position of the wavefront 
is the zero level set of <p at time t as defined by 

4 +F|V0(x,Ox| =0. (2) 
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[033] If the speed function F does not depend on time and preserves its sign, 

then the level set function <p(x, 0 = 0 becomes single-valued in time t, i.e., each 

pixel is evaluated only one time. This leads to the well-known Eikonal equation 

|V^(x)| = F, (3) 

where the hyperfunction now corresponds to an arrival time of the wavefront at the 

point x. One of the subtleties of this equation is that the solution may be non- 

differentiable. 

[034] In 91 2 , a finite difference approximation of the gradient term in the 

Eikonal equation (3) can be obtained as 

| V yr(x) | 2 « max( D*yr , 0) 2 + min( D'yr , 0) 2 + 

max( /);> , 0) 2 + min( D'y/ , 0) 2 = F 2 , (4) 

where the difference is selected in terms of the direction of the flow information as 

U x Y — A *> U x V A ' \?) 

Ax Ax 



[035] Further simplification can be achieved by using a uniform grid of pixel 
points, i.e., Ax= Ay=l. The above formulation is called an 'upwind scheme'. It is 
straightforward to prove that this scheme converges to the correct viscosity 
solution. The upwind scheme means that information propagates one way, from 
smaller values, to larger values of y/ . 



[036] Fast Marching 



[037] A fast marching method solves the Eikonal equation, which applies in the 
case of a convex, non-negative speed function. Because the wavefront crosses each 
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grid point only once, the moving wavefront problem can be converted into a 
stationary formulation in which the arriving wavefront contains the information 
about what is moving. 

[038] This can be expressed by Huygens Principle: "Every point on a primary 
wavefront serves as the source of spherical secondary wavelets such that a primary 
wavefront, at some later time, is the envelope of these wavelets." 

[039] To give an example, consider building a brick wall. One places a brick on 
a current level of the wall, and then moves laterally, at the same level, and puts 
another brick on the same level. After all the bricks are placed at a given level, one 
goes up to the next level and repeats the process. The important thing is that the 
wall is built from the ground up, each level must be completed before the next 
level is started. This is, in fact, the main difference between the prior art linkage- 
based region growing method, and the level set method according to the invention. 
The linkage-based region growing does not check whether a level is completed 
before moving to the next level. The linkage-based method can add bricks to the 
wall at any place where there is an empty space around a brick. This can surely 
cause the wall to extend in an uncontrollable manner and tumble down. 

[040] In contrast, the fast marching method, as described above, imitates the 
wall building process. Starting with an initial position for the wavefront, the 
method systematically marches the wavefront outwards one grid point at a time, 
relying on entropy-satisfying schemes to produce the correct viscosity solution. 
Similar to the brick-laying process, an Eulerian algorithm determines the arrival 
times of the wavefront starting from an initial wavefront 0. 
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[041] Alive, Narrow Band, and Far Away Sets 

[042] An alive set includes all image points (pixels) y/ {x,} such that the arrival 
time of a point ^"(x,) = 0 of the wavefront is zero. Initially, the alive set 
corresponds to {x, : y/ (x,) = 0}. 

[043] A narrow band set includes neighboring image points (pixels) {x ; } of 
points in the alive set. The arrival time of the points in the narrow band set is y/ 

(xj) =1= F(xj). In essence, the narrow band set of points defines the wavefront 
according to the invention. 

[044] The remaining points in the image are in afar away set of points. 

[045] Then, the wavefront is moved forwards using the following iterative 
process. 

Let x* be a point in the narrow band set with a smallest arrival time 
^(x*); 

Move point x* from the narrow band set to the alive set; 

Move all neighboring points x j of the point x from the far away set 
to the narrow band set if they are not in the narrow band set; 

Update y/( x* 7 ) for all neighbors of x*, i.e., solve the Eikonal equation 
using the quadratic approximation given in the equation (4); and 

Repeat until no point remains in the far away set. 
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[046] A heap-sort technique can be used to locate systematically the proper 
grid point to update. Thus, it is unnecessary to backtrack over previously evaluated 
points. The resulting technique processes N points in N logN steps. 

[047] Segmentation 

[048] As shown in Figure 1, the above wave propagation method is used to 
partition 100 an image 101 into homogeneous segments 109. Base points 111 are 
selected 1 10 in the image 101. Wavefronts 121 are initialized and emitted 120 
from the selected base points 111. The wavefronts 121 are propagated 130 
according to a speed function. Then, properties of these wavefronts are evaluated 
to determine 140 if a termination condition 141 is satisfied. If the termination 
condition is satisfied, then the final position of wavefronts or 'boundaries' are used 
to produce 150 the segments 109, otherwise the propagation continues. 

[049] In other words, each base point 111 serves as a seed of the segment to 
which the base point belongs, and the wave propagation 130 expands the segment 
by wavefront propagation until the termination condition is satisfied. 

[050] As shown in Figure 2, further details of the method according to the 
invention can be summarized generally as: 

Select 110 a base point 111; 

Initialize and emit 120 a wavefront 121; 

Propagate 130 the wavefront 121 according to a speed function 131; 
Determine 140 if a termination condition 141 is satisfied, see Figure 1; and 
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Mark 135 the current segment, and remove the current segment 109 from set 
of available points. If no point remains in step 140, then stop, otherwise iterate 
with step 1 to select the next base point. 

[051] In Figure 2, the termination condition is a check to determine if any pixel 
that is not included in any of the previously segmented regions remains in the 
image. If there are such pixels, then another base pixel is selected 110, and the 
wavefront propagation 130 is repeated. 

[052] In other words, the step 140 in Figure 2 represents a process that is 
applied pixels are marked due to wavefront propagation result, which leads to a 
segmented region. The marking 135 applies a distinct label to the pixels that are in 
the same region that is obtained by wavefront propagation, i.e., all pixels in the 
region have the identical label, e.g., a number, which is different from the labels of 
other regions. 

[053] Selecting Base Points 

[054] The base points 1 1 1 are selected 110 iteratively, that is one at the time. 
After a base point has been selected, the wavefront is propagated, and the segment 
is identified. An image point (pixel) is selected as a base point b if the point 
represents a distribution, e.g., a color or intensity distribution, within a relevant 
local neighborhood. A variance a , and a gradient V/ are good indicators of color 
or intensity homogeneity. Therefore, construct a gradient image 102 and a variance 
image 103 in preliminary step. 
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[055] The likelihood of being a base point is inversely proportional to the 
gradient and variance. A point has a high likelihood of being a base point if it a 
least gradient value and a least variance value with respect to global maximum 
gradient and variance values. The gradient and variance indicate that the color of 
the point is consistent with neighboring points, and that the color distribution is 
uniform around the point. This selection process prevents initializing and emitting 
a wavefront at an edge, which can disturb the local homogeneity of the segment. 
The selection 1 10, as described above, effectively centers the initial wavefront for 
a more accurate segment boundary localization. Ideally, a base point should be at 
the center of the segment to which the base point belongs. 

[056] Let I(x, c) be the color value of the point in a corresponding color 
channel c. Assuming the color space is orthogonal, e.g., RGB, the local color 
variance within a local window W around point is defined as 



where ju (x, c) is the mean of the color within the window, and co is an envelope 
function, which has a Gaussian form. A contribution of all points within the 
window is same, i.e. co = 1. The color gradient is 



[057] Using the gradient and variance images 102-103, a dissimilarity score 
d(x) is determined for an image point as 



( x > = I 2>['(x,c) - Mx,c)] 2 , 



(6) 



| V/(x) = XI /(*) + 1, y, c) - I(x - 1, y, c)\ + 



\l(x,y+\,c)-I(x,y-\,c)\. 



(7) 




(8) 



k 
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where K - 1 is the number of the previously emitted wavefronts, k is the arrival 
time of the k th wavefront, and#, , a 2 , and a 3 are constant weights. 

[058] Initially, there is only one wavefront K = 1 , and the arrival times y/ (x) 
are set to a maximum number. For the Eularian method, the contribution of the 
previous arrival times is dismissed, i.e., the weight is set to zero, i.e., or, = 0, and 
only the points that are not included in the previous segments are considered. For 
the non-Eularian method, all image points are evaluated. The point with the least 
dissimilarity score is assigned as the next base point according to 
b K : argmax<i(x). (9) 

[059] Figures 4-6 show sample images 401, 501 and 601, and their 
corresponding local color variance images 402, 502, and 602. To determine the 
variance within a bigger window without increasing the computational load, 
equation (6) is applied using the same window size (5 X 5) in a hierarchical manner 
at different image resolutions, e.g., the original resolution, the original resolution 
downsampled by two, and the original image downsampled by four. 

[060] The mean and variance images are constructed at each resolution level. 
Then, the mean and variance at each level are aggregated by compensating for the 
downsizing. For each base point, a substantially circular wavefront K is initialized 
and centered on the base point. The wavefront is 'substantially' circular in the 
sense that it is understood that pixels in an image are discrete and regularly spaced. 
Initially, only the base point bK and its 4-neighbors are on the wavefront K. The 
points on the current wavefront also correspond to the narrow band B K set defined 
above. 
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[061] A mean color Sk of a segment is assigned to the current segment K;S K = 
I(b K ). 

[062] Wave Speed 

[063] The conventional Eikonal equation assumes that the speed of the wave on 
the normal direction to the wavefront is constant. In contrast, in the method 
according to the invention, the speed of the wavefront is modified as a function of 
color to make the propagation adaptive to the underlying color distribution. 

[064] Therefore, three different sources of information are embedded into the 
speed function, namely a color similarity of a current point to a current segment, a 
color gradient, and a color variance. In case the color value of the current point in 
the narrow band is similar to the segment's color mean, the speed is increased to 
include the point in the segment. If a point has significantly different color value, 
the speed of the wavefront propagation decreases at this point. The speed function 
is set as a function of a color distance 
F(x) = g||I(x), S K ||, |/(x)|, o- 2 (x), (10) 
where ||.|| is a distance norm, e.g., a magnitude of the difference (Li). 

[065] A simple speed function is given by F(x) = |/(x) - S K \. At each iteration of 
the fast marching, update the mean S K of the segment color is updated, the color of 
the narrow band set B K , the number of points N K in the alive set of points, and the 
number of points M K in the current narrow band set as 
S' K =l/N t K [N' l K S K + I(x*)], 
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&K = \IM tK [M' \B K - I(xj) + (xj )], 



N t K =N ,1 K +l,and 
M K = M"- l K -1 + C 



(11). 



[066] The speed of the wave is determined by using the above color mean 
parameters to achieve color homogeneity of the segment. 

[067] Termination Conditions 

[068] An important problem of region growing is determining termination 
conditions. Four different termination conditions are possible for the Eulerian 
approach, as shown in Figure 3. The conditional expression for a particular 
termination condition is determined in step 124. One condition measures a color 
similarity 301 of the segment and the current point according to 



[069] That is, if the color value of the current point is greater than that a color 
similarity threshold t a 141, then the wavefront stops at that point. The threshold 
141 can be set to a constant. A threshold r A = 3.2 m /64 produces accurate 
segmentation for most images under the assumption that color channels are m-bit 
coded. This value can be obtained empirically. 

[070] Another measure makes a comparison between the point and the 
segment, and also compares the whole wavefront with the segment according to 



Here, the threshold is similar to the previous color threshold, i.e., t a = t b . 



\\S K -I(x)\\<T A . 



(12) 



\\S k -B k \\<t b 



(13) 
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[071] The segment can also be expanded up to the edges: 
|V/(x*)<r c (14) 
However, this termination condition 142 relies on the assumption that the speed of 
the wavefront is constant. Thus, the wavefront can continue moving until it reaches 
an edge. In this formulation, a point that has a larger gradient value than the 
threshold T c 142 is considered as an edge point. In case the wavefront reaches an 
edge, the wavefront stops marching at that point, and the propagation continues 
until no more expansion is possible. Thus, this condition depends on a gradient 
similarity 302. 

[072] Alternatively, for a non-constant wavefront speed function, the arri\^l 
time y/(x) of the wavefront can used as a termination condition 303. The wave 
continues propagating as long as \ff{x) <T D . (15) 

[073] It is necessary to consider the effect of the size of the image size and the 
parameterization of the speed function when the threshold is assigned. If the arrival 
time threshold z D 143 is not scaled with the image dimensions, a wavefront can 
stop prematurely in case of a large segment. 

[074] In Figure 3, the termination step 140 checks a conditional produced by 
step 124 to determine whether to continue wavefront propagation 122, or to 
terminate using one of the similarity scores determined in steps 301-303 and one of 
the corresponding thresholds 141-143, e.g., color similarity with color threshold, 
gradient magnitude with gradient threshold, time with arrival time threshold. The 
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termination condition is score > threshold. If the termination condition is not 
satisfied, then the wavefront is updated 125, and the parameters are updated 126. 

[075] Thus, the termination condition in step 140 of Figure 3 is a conditional to 
determine whether to propagate the current wavefront, while in Figure 2, the 
termination condition determines whether to initialize a new wavefront. 

[076] Speed Function 

[077] The method according to the invention can be used with a constant or 
varying speed function 

F(x) = a x \s k - 7(x)| + a 2 | V/(x) +a 3 , a x > a 2 > « 3 . The constant term a 3 
incorporates the effect of the distance from the base point. Thus, the term serves as 
a reinforcement on the shape of the wavefront. In other words, the constant term 
acts as a viscosity parameter and controls the curvature of the wavefront. Note that, 
in case of a 2 » a 3 , the wavefront moves forward at a constant speed by 
disregarding the color similarity of the wavefront to segment. 

[078] Figure 4 shows the motion of the wavefronts 401 and the corresponding 
arrival times 402 for the image of Figure 2. As is visible from the Figure, the 
wavefront can successfully modify its shape to match the underlying color 
distribution. 

[079] Figures 5-7 show the extracted segments for the corresponding images in 
Figures 1-3 using the method according to the invention. 
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Effect of the Invention 

[080] The invention provides a method for partitioning an image into segments 
by propagating a wavefront. The invention overcomes limitations of prior art 
region growing techniques. The invention uses a color similarity measure, a 
gradient based speed function, and four different propagation termination 
conditions. 

[081] Using the level set techniques, the invention effects concurrently 
integration of spatial continuity, boundary shape control, and color consistency. 
This is an important improvement over conventional linkage-based region growing 
methods. 

[082] Furthermore, the fast marching solution of the Eikonal equation 
effectively reduces the computational cost of the segmentation process. The 
method is robust and accurate. 

[083] It is to be understood that various other adaptations and modifications 
can be made within the spirit and scope of the invention. Therefore, it is the object 
of the appended claims to cover all such variations and modifications as come 
within the true spirit and scope of the invention. 
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